Distribution of supercurrent switching in graphene under proximity effect 
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We study the stochastic nature of switching current in hysteretic current-voltage characteristics 
of superconductor-graphene-superconductor (SGS) junctions. We find that the dispersion of the 
switching current distribution scales with temperature as aj oc T'^'^ with aa as low as 1/3. This 
observation is in sharp contrast with the known Josephson junction behavior where aj oc ' with 
a J = 2/3. We propose an explanation using a generaUzed version of Kurkijarvi's theory for the 
flux stabiHty in rf-SQUlD and attribute this anomalous effect to the temperature dependence of the 
critical current which persists down to low temperatures. 



PACS numbers: 74.45. +c, 72.80.Vp, 74.40.-n, 74.50. +r 

Since the extraction of single-layer graphene [H [5] 
much effort has concentrated on its study due to the 
promising potential in applications. The knowledge of 
graphene properties and expertise in making high qual- 
ity devices have grown substantially [SHS]. Neverthe- 
less, the transport in graphene subject to nonequilib- 
rium conditions and in the proximity to a superconduc- 
tor, an important ingredient in the majority of applica- 
tions, is far from being fully understood. Unlike metal- 
superconductor interfaces reflection from a graphene- 
superconductor boundary is governed by the specular 
Andreev processes [Sj . This peculiar effect combined with 
the unique band structure of graphene makes proximity 
effect in graphene a particularly interesting subject to 
study. 

Recent experiments on the superconductor-graphene- 
superconductor (SGS) devices have revealed many inter- 
esting features caused by the proximity effect 1^. These 
include an observation of supercurrent and subsequent 
measurement of the current-phase relation, signatures 
of multiple Andreev reflection in the differential con- 
ductance, and Shapiro steps under microwave irradia- 
tion, see Refs. [8til2j. Recent measurements have also 
revealed the residual resistance of SGS junctions for cur- 
rents below critical which was attributed to the phase 
diffusion phenomenon ^13] followed by the crossover to 
macroscopic quantum tunneling regime at low tempera- 
tures 14 . Here we report first systematic study of ther- 
mally activated dynamics of phase slips in SGS junctions 
through the measurement of the switching current distri- 
bution. 

Measurement of the decay statistics of metastable 
states is a powerful tool for revealing the intrinsic ther- 
mal and quantum fluctuations. In the Josephson junc- 
tions ( JJ) a metastable dissipationless (superconducting) 
state decays into dissipative (phase slippage) state when 
the bias current / reaches a critical value called switching 
current Iswj which is stochastic. Analysis of the distri- 



bution of the switching current was employed to reveal 
macroscopic quantum tunneling in JJs jl5j . supercon- 
ducting nanowires [TS] , small underdamped JJs [T7] and 
intrinsic JJs in high-Tc compounds |18j . Experimentally 
observed temperature dependence of the switching cur- 
rent dispersion aj always follows a power law aj oc T^/^ , 
if the switching is induced by a single thermally acti- 
vated phase slip (TH]. However at sufSciently low tem- 
peratures the temperature dependence of ui saturates, 
which is usually attributed to the macroscopic quantum 
tunneling |15| . 

Study of switching current distribution in conventional 
SNS junctions, where N is normal metal, is obstructed by 
the fact that such junctions are usually overdamped. As a 
result their I-V characteristics are smooth and the notion 
of the switching current is not applicable. Here we report 
a study of moderately underdamped SGS junctions with 
the quality factor Q ~ 4 for the entire span of gate volt- 
ages. Our main finding is the anomalous temperature 
dependence of the switching current dispersion cr/ cx T"'^ 
in SGS devices with 0.3 < aa ^ 0.5, which persists for 
a wide range of gate-induced doping and is significantly 
smaller than the usual aj = 2/3. In general, any power 
law different from 2/3 is associated with the possibility 
of quantum phase slips. In our graphene-based proximity 
junctions, although the power law notably deviates from 
2/3, we argue that thermally activated phase slips are the 
major contributor. We interpret an anomalous dispersion 
of a I by using a generalized Kurkijarvi model ^20j . Our 
conclusion is that the slowed temperature scaling of c/ 
in SGS junction is due to the substantial temperature 
dependence of the critical current, which persists down 
to low temperatures in SGS systems. 

Graphene flakes are deposited on 280 nm thick Si02 
surface using mechanical exfoliation |lj. Raman spec- 
troscopy is used to confirm the number of layers |21j . 
The electrodes, which have a fingered shape (Fig. [l|i), 
are patterned from a bilayer Pd/Pb (4nm/100nm), as 
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FIG. 1: [Color online] a) SEM micrograph of sample 105. 
Distance between the electrodes along the current (length of 
the junction) is L = 265 nm. Width of the junction (distance 
across the current) is W — 214 /im. For sample Ills, L = 
280 nm and W — 9.9 fim. b) The hysteretic I-V curves of 
SGS junction (sample 105) taken at various gate voltages. 
The switching Isw and retrapping Irt currents are shown. 

explained in the supplementary materials (SM). In order 
to measure the switching current distribution, the ampli- 
tude of the sinusoidal current bias is set somewhat higher 
than the maximum switching current, and it is adjusted 
when needed to keep the sweep speed roughly constant. 
The number of switching events for each distribution was 
either 5000 (for the sample Ills) or 10000 (for the sample 
105). At low temperatures the I-V curves of the samples 
exhibit a hysteretic behavior. Fig. [TJa, which enables us 
to study switching current statistics. 

Our main focus is on the cr/(T) function. Figure [2] 
shows a switching histogram for sample 105 at (a) Dirac 
point {Vg = —30V) and (b) Vg = 50V. During the exper- 
iment, some anomalously premature switching events are 
recorded. These events, which significantly deviate from 
the general population of the distribution, are very rare 
and are believed to be unrelated to thermal fluctuations. 
In order to exclude these anomalous jumps from the stan- 
dard deviation calculation we first convert the raw data 
to the switching rate r(7) according to the Kurkijarvi 
method [121 HO] ■ The Kramers and Stewart-McCumber 
theories combined (see below) lead to the expectation 
that InT cx (l-///c)^/^ In Figs. [2]: and[2|l we plot InF 
versus (1 — I /IcY^'^. The critical current Ic is tuned 
to make the graph as linear as possible. Then the lin- 
ear part of the graph is fit with a straight line. Hollow 
squares and circles are the measured data. Filled sym- 
bols are those points which were used to find the best 
linear fits. The best fit r(/) is then used to regener- 
ate the distribution of Isw by inverting the Kurkijarvi 
transformation. The results are shown as black curves 
in Fig. |2^ and Fig. [2]d, for the experimental sweep rates 
were 363/xA/s and 2.7mA/s. The red curves are com- 
puted distributions for dl/dt = 100/zA/s and l.OmA/s 
correspondingly. The dispersion cr/ is then computed for 
the same value of dl/dt for all temperatures. 

Our main results are presented in Fig. |3] This fig- 
ure shows standard deviation aj and critical current Ic 
versus temperature for various gate voltages. Figures |3^ 
and[3]D are log-log plots of cr/ versus T. The best linear fit 



provides aa, which is defined by the equation aj oc T°"^ . 
The estimated error or uncertainty in the power values 
is about 7%. Overall, the best fit ac's are different from 
the theoretically predicted JJ value aj = 2/3 — 0.667. 
Since numerous previous experiments on JJs established 
the power close to 2/3 while our data indicate powers 
roughly between 1/3 and 1/2, an understanding of such 
discrepancy is desirable. 

We interpret these observations based on the follow- 
ing model. Since the pioneering theoretical work of 
Kurkijarvi [20^ and its experimental confirmation by Ful- 
ton and Dunkleberger [19] kinetics of stochastic phase 
slips in the JJs is described within Stewart-McCumber 
model 23J, which employs sinusoidal current-phase rela- 
tion (CPR), /s(0) = Icsm{(j)), and represents the total 
current as a sum of superconducting, normal and dis- 
placement components. At the mesoscopic scale and, in 
particular, in the context of graphene proximity circuits, 
there are reasons to question the applicability of such 
model given the possibility of a highly nontrivial struc- 
ture of /s(0) (see SM). This naturally raises a question 
about the universality of the previous results with respect 
to the form of the CPR. It is rather remarkable to realize 
that the predictions of the theory j20j in fact extend be- 




FIG. 2: [Color online] (a), (b) Switching current distributions 
at Dirac point (Vg — —30 V) and Vg = 50 V. The black curve 
shows a theoretical fit to the experimental distribution with 
an experimental speed: (a) 363/iA/sec and (b) 2.7 mA/sec. 
The red curve shows a calculated distribution with a new 
sweeping speed: (a) 100/iA/sec and (b) 1.0 mA/sec, using 
the escape rate to the standardized sweep speed, (c)-(d) 
Logarithm of the escape rate is shown as a function of the 
scaled current. The raw data are shown as hollow circles and 
squares. Filled squares and circles are used to calculate the 
critical current and to fit the escape rate, which is shown as a 
solid line. Anomalous premature jumps are visible as isolated 
data on the left side of the graph. 
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FIG. 3: [Color online] For sample Ills {Vgo = -IV) data 
sets on panels (a), (c) and (e) correspond to gate voltages 



to Vg=50, 30, 10 



30 V) data sets on panels (b), (d), (f) correspond 
10, -30, and -50 V, from top to bottom. 
(a)-(b) Standard deviation vs. temperature, in the log-log for- 
mat. The best linear fits determine the power ac, which is 
shown near each fit. (c)-(d) Critical current vs. temperature 
for sample 105 and Ills at various gate voltages. Solid lines 
are theoretical fits [55]. (e)-(f) Normalized standard devia- 

1/3 

tion, oi/Iq oc T°"-^ , vs. temperature, in the log-log format. 
The corresponding powers aa are indicated. The graphs are 
shifted vertically for clarity, (g) Log-log plot of standard de- 
viation vs. critical current at four different temperatures, (h) 
Log-log plot of the scaled standard deviations vs. scaled crit- 
ical current. 

yond the limits of its original validity. We now proceed to 
the generalization of the Kurkijarvi's theory [2D] devel- 
oped for the statistics of thermally activated phase slips 
in a flux-biased rf-SQUID to the case of a current-biased 
weak link with an arbitrary CPR. 

Within the Stewart-McCumber model the dynamics 
of the phase 4> is equivalent to the dynamics of a vis- 



cous Brownian particle subject to the following external 
potential: 



G(0) = - M0/2e, 



(1) 



which is the Gibbs potential. Here F{(f)) — 
{h/2e) J d(j)Is{(j)) is the free energy and / is the bias cur- 
rent. We assume that Is — Isi<P) is a single- valued 
smooth function. For / = 0, G{cf)) is a periodic func- 
tion of (j) with alternating local maxima and minima. In 
the absence of fluctuations the phase is trapped in one of 
the minima as long as I < Ic, which is a state with zero 
voltage. In the resistive state when I > Ic the phase in- 
creases with time. In the presence of thermal fluctuations 
even at I < Ic the phase can escape its local minimum, 
i.e. experience a phase slip, which drives the junction 
into a phase-running resistive state. Such a process is 
detected as a voltage jump (a switching event) on the 
I-V curve. Upon decreasing / the junction may show 
hysteretic behavior and the quality factor Q determines 
the width of the hysteretic region. The activation rate 
of a phase slip for a moderately underdamped {Q > 1) 
to overdamped {Q <C 1) junction is given by Kramers 
theory [24: (hereafter fc^ = 1) 



r = (l/2^)(VW4T^-?7/2)exp(-AG/r). (2) 

The energy barrier AG is the spacing between two con- 
secutive extrema: AG = G{4>+) — G{(f)-). The prefactor 
is determined by the curvature of the potential at mini- 
mum, = C~^(2e/h)^d^G and by the damping param- 
eter r] = l/RffC, where Rn and G are effective normal 
resistance and capacitance of the junction. Notice that 
Q = LOp/r] where uip = •\/2e/c//iG is plasma frequency. 

To find the activation barrier AG(/) let us introduce 
a critical phase <pc defined through Ic — Is{4'c)- In 
the vicinity of (pc one can use Taylor expansion Is {(p) — 
Ic — \\I'l;\{<p — cpcY provided /s{(/)) is a smooth func- 
tion. F{(l)) is obtained by integrating the supercurrent 
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^ ^'^I{^)\<Pc- These equations determine the loca- 
tions of the two consecutive extrema of the Gibbs poten- 
tial: (f>±-(j)c = ±^2{Ic- I)/\I'S\. Using Eq. Q one 
finds 



over the phase, which gives F{(j)) = Fc + ^Ic{4' — 
^c) - ^^l^cK-^ - '^c)^ where Fc = F(0c) and 



AG(/) = Gc(l-///c)'/', Gc = 




(3) 



The curvature of the Gibbs potential at the extrema 
points can be obtained in a similar way and is given by 

u:'^^^p^2i\I'S\/Ic)il-I/Ic)- 

The knowledge of the decay rate (Eq. [2| allows one to 
determine the probability p that a phase slip occurred by 
the time t, which reads: p{t) = 1 — e^-^ '"f* , where 
F = r{I{t)). Note that the probability of not having a 
phase slip by the time t is 1 — p. For a constant bias 
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current sweep Wi = I^^dl{t)/dt the probability p can be 
evaluated analytically. Introducing reduced current vari- 
able u = [Gc/TY^^{1-I /Ic) and recalling the definition 
of the quality factor one obtains the following expression: 

_j,^-.3/2 2Tujll{iTTu:,T^Gc) 
p=l-e , X = 

1 + v/l + Q'(r/Gc)l/3ul/2 

(4) 

This result for the probability of a phase slip holds for a 
moderate to high damping provided Gc ^ T , the condi- 
tion which is very well satisfied in most of our measure- 
ments. In the limit of high damping when Q — > and 
X — ?> Tujp/{3TTUJiriGc) one recovers the result of |20) . 

To evaluate the dispersion of the switching current we 
notice that the probability distribution P{x) of a vari- 
able X is obtained from p{x) by the differentiation with 
respect to x i.e. P{x) = —dp{x)/dx. Given the relation 
between the bias current / and the reduced current u 
this implies that the dispersions of these variables are re- 
lated as (7/ = Ic{T /GcY^'^au- The crucial observation is 
that dispersion cr„ considered as a function of X is con- 
stant within a few percent while X is varied by several 
orders of magnitude, so that for all practical purposes 
Gu is temperature independent j20j. Using Eq. ([s]) and 
assuming that the temperature scalings of Ic and I'^ are 
the same we obtain the following temperature scaling for 
the dispersion of switching current: 

ai{T)^{T/'fof'Hy\T), (5) 

where $0 = h/2e is the flux quantum. Eq. ([s]) is the main 
result of the calculation, which describes the tempera- 
ture dependence of aj for any smooth CPR. According 
to Eq. ([5]) the power of 2/3 in the temperature scahng for 
a J is only expected if Ic{T) = const. In the SGS junc- 
tions the critical current keeps increasing down to the 
very low temperatures (Fig. |3j: and [sji) , due to the diver- 
gence of the normal metal coherence length in graphene, 
thus leading to the stronger proximity effect. The solid 
lines in the figures are the fits to the SNS junction the- 
ory [2S]. The fitting parameters for the theoretical fits 
are mean free path /e=10-25nm, which is similar to pre- 
viously reported values [MTU] , and the normal resistance 
Rjq , which is of the same order of magnitude as the one 
measured directly. 
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To further confirm our conclusions we plot aijl^J, 
versus the temperature as suggested by Eq. ([s]). The 
results are shown in Fig. and |3f. Such critical- 
current-normalized dispersion obeys the power law with 
the power close to 0.6. For sample 105, the power rests 
between 0.47 and 0.62; for sample Ills these values vary 
between 0.55 and 0.63, which are close to 2/3, predicted 
by the adopted Kurkijarvi model. 

Another type of scaling which is suggested by Eq. ([5| 
and which can be accessed experimentally is the depen- 
dence of the dispersion on the critical current at con- 
stant temperature. Doping-dependent conductivity of 



graphene provides a unique possibility to vary the criti- 
cal current while keeping the temperature constant - an 
experimental "knob" which is inaccessible for other types 
of junctions. 

In Figs. [3|; and[3)i we present results of such measure- 
ments. We plot the ln(CT/) versus In(Jc) for various tem- 
peratures (Fig. [3^) and scaled ln(cr7) versus scaled ln(/c) 
(Fig.jsjr). The average value of the power for sample Ills 
is 0.34. The power of scaled data for sample 105 is 0.38, 
which is shown as the best fit of the data. The resulting 
powers are very close to the theoretically expected value 
of 1/3. 

In summary, we have studied the dispersion of the 
switching current distribution in moderately under- 
damped SGS junctions with clear hysteretic I-V char- 
acteristics. A systematic measurements of both temper- 
ature and critical current scaling (at constant T) of the 
dispersion is performed. The latter study, unavailable 
in regular junctions, is made possible by a gate-voltage- 
tuned conductivity of graphene. The temperature scaling 
of the switching dispersion shows unusual power laws, 
which is explained theoretically by taking into account 
the temperature variation of the critical current. The 
critical current scaling of the dispersion is explained theo- 
retically by combined Stewart-McGumber and Kurkijarvi 
models, and is applicable for the mesoscopic junctions 
with arbitrary current-phase relationships. 
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Supplementary Materials 

Raman spectra and additional sample preparation 
details.- In Fig. [4j the Raman spectroscopy results are 
plotted for two different thin graphite flakes. In the 
thicker flake, the peak profile is non-Lorentzian and 
shows a small side peak. In addition, the FWHM of 
the best fit Lorenzian is 33 cm^^. However, in the thin 
fiake, the peak profile is Lorentzian and located between 
2600 and 2700 cxnT^ . In addition, the FWHM for this 
sample is 13 cm~^, which is suggestive that the thin flake 
is graphene [T]. 




2550 260(1 2650 J 2700 2550 2600 2650 _j 2700 

Raman Shift (cm ) Raman Shift (cm ) 



FIG. 4: [Color online] (a) Raman spectroscopy for a thick 
graphite flake showing a double peak. Inset: optical image 
of the thick graphite flake. The red dot shows where the 
data are recorded, (b) Raman spectroscopy of a thin graphite 
flake showing a single peak. The peak is fit to a Lorentzian 
curve centered about 2628.5 cm"^ with FWHM of 13 cm~^ 
suggesting that it is graphene. Inset: optical image of the 
graphene flake. 

A pair of pseudo-four-probe electrodes are patterned 
using ebeam lithography. In a thermal evaporator, a 
layer of 4nm Pd is evaporated at a rate of 0.5-1.0 A/s 
as slower evaporation rates would lead to a high contact 
resistance, possibly due to heating of the graphene by 
radiation coming from the molten target. Following the 
first layer, a layer of 100 nm Pb is evaporated at a rate 
of 10-30 A/s. We have found that a faster Pb evapo- 
ration provides more uniform films, thus it is desirable. 
We have used a liquid nitrogen trap to remove the resid- 
ual contaminants in the chamber and thus to increase 
the quality of the films. After evaporation, the sample 
is placed in a hot acetone bath for total 5 minutes for 
the lift-off step. While being in the acetone bath, the 
sample is sonicated for 10 seconds every other minute. 
The sample is kept in the second bath of acetone for 5 
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minutes before it is rinsed with iso-propanol and dried. 
The Pb electrodes are too soft to use wire bonder thus 
indium dots have been used to connect the leads to the 
chip. The sample is measured in a Hea system. The sys- 
tem is equipped with room temperature of Pi Filters and 
low temperature Cu powder filter. 

Remarks on junction capacitance.- Total capacitance 
C is an important characteristic of the junctions. It 
defines quality factor Q which in a way determines 
whether junction's dynamics is overdamped or under- 
damped. In the main text we determined Q from the 
experimental values of the critical and retrapping cur- 
rents, which together with the junction's normal resis- 
tance provided us an estimate for C. An alternative 
way to determine C is suggested by considering the 
geometry of our sample, see Fig. [5] An estimate of 
the capacitance across the graphene sheet can be found 
with the help of the formula for coplanar electrodes [2J: 
Cg — ereoWK{k)/K{\/l — k"^), which leads us to the 
conclusion that the electrode capacitance should be ^ 10 
fF. Here eo is the permittivity of free space, er is the 
relative permittivity of silicon ~ 4, is the junction 
width ~ 300/im, and K{k) is the complete elliptical 
integral of modulus k = (1 + L/W)~^; L is the dis- 
tance between electrodes. Note, however, that two elec- 
trodes have an additional capacitive coupling through 




FIG. 5: [Color online] (a) The capacitance across the 
graphene, Cg, can be calculated using the cross sectional area 
of the Pb film and the distance between the Pb electrodes, 
(b) The Pb electrodes form capacitors through the Si02 with 
a plate separation of d. (c) The electrical scheme of the capac- 
itance through the graphene and the two capacitors formed 
between the electrodes and the gate. The equivalent capaci- 
tance is Cg + CelCe2/(Cel -|- Ce2). 



6 




current as the function of applied magnetic flux can be 
fit to a Fraunhofer function: 



-10 10 
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FIG. 6: [Color online] A color plot of the differential resis- 
tance versus bias current and magnetic field. The red color is 
zero differential resistance while the green is the normal state 
resistance. This plot also shows the expected Fraunhofer pat- 
tern. 

the gate. Thus Fig. [5]d gives the proper geometry and 
electrical scheme for the calculation of the full capaci- 
tance across the junction. In this geometry, d is 280nm, 
and the area of each electrode (including the presence 
of a pressed indium dot to facilitate the electrical con- 
nection to the pins on the chip carrier) is 2.9 and 2.27 
fim^. Thus the capacitance of each electrode to the gate 
is 36.7 and 28.7 pF, which is calculated from the formula 
Ce — Cr^oA/d, where d is the distance to the substrate. 
Using the circuit scheme as in Fig. [5];, the total capaci- 
tance C = Cg+ CelCe2/(Cel + Ce2) IS fouud tO bc 16pF, 

which is in between the values that can be determined 
from the quality factor argument. Indeed, for 330 mK, 
the quality factor was found to vary between 3.5 and 4.2 
for the entire range of gate voltages from the values of Ic 
and InT- Thus the SGS junction (sample 105) is mod- 
erately underdamped as it is larger than 0.84. Since the 
quality factor is Q 



WpRnC where Wp — y/2eIc/fiC, 



we calculate that the capacitance of the sample varies 
between 12 pF and 50 pF. This consideration proves 
consistency of our estimates between two independent 
approaches, and confirms the conclusion that our SGS 
junction is in the underdamped regime. 

Characterization of graphene proximity junctions - In 
this section, the basic characterization of the electronic 
transport properties of graphene proximity junctions will 
be presented. We focus mainly on the magnetic field de- 
pendence of the supercurrent, which demonstrates high 
quality JJ in our graphene based devise, and on the gate 
voltage dependence of resistance and switching current, 
which traces properties of the graphene, associated with 
the Dirac spectrum. This discussion complements pre- 
sentation given in the main text. 

The magnetic field dependence of the switching current 
can be measured and compared to Josephson junction 
theory. In an extended J J, the maximum (or critical) 



m 

/(O) 



sin(7r$/<I>o) 



7r<i)/$c 



(6) 



where $ is the magnetic flux through the junction which 
is given by $ = BW{L + A/2), B is the magnetic field, 
and $0 — h/2e is the magnetic flux quantum. The ef- 
fective area is larger than just LW due to the magnetic 
field penetration into the electrodes by a distance equal 
to the penetration depth A. 

The fit of our data to the Fraunhofer pattern can be 
seen in Fig. [T]l of the main text, where a good fit is ob- 
tained using an area that is estimated from the junction 
area of ^ 10~^° m^, which results in a magnetic field 
period of ~ 20 mT. A vertical shift of 2.5 mA was also 
used to obtain the fit to account for the supercurrent 
observed even at the Dirac point. The magnetic field 
period obtained from the experimental data is somewhat 
small at about 7 mT. This discrepancy between the fitted 
and measured value can be traced to an understanding 
of the effective area of the junction W{L + A/2). Using a 
penetration depth of 120 nm results in an effect area of 
~ 2 X 10~^° m^ and a magnetic field period of 9 mT. If the 
area of the entire electrode plus junction were used, the 
magnetic field period would be ~ 7.8 mT. Thus, to fit to 
the Fraunhofer period the penetration depth in our sam- 
ples should be larger than 120 nm in order to increase the 
effective area of the junction. The differential resistance 
of the junction can also be measured as a function of bias 
current and magnetic field, see Fig. [6j Again, the Fraun- 
hofer pattern is observed, with a slightly reduced value 
of the switching current even at a applied gate voltage of 
50 V due to a shift of the Dirac point after the thermal 
cycling. 

The switching current and resistance as a function of 




-50 -40 -30 -20 -10 
V„ 



10 20 30 40 50 



FIG. 7: [Color online] (a) The normal resistance, the retrap- 
ping current, and the switching current shown as a function 
of gate voltage. As the gate voltage is swept from the Dirac 
point, the switching current and retrapping current increase 
however resistance of the sample decreases. 
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the gate voltage are displayed in Fig.jTj As the gate volt- 
age is swept away from the Dirac point, the total carrier 
concentration increases and thus, the switching current is 
expected to increase as well, which is confirmed in Fig. [7^. 
Similarly, the retrapping current increases as well. In ad- 
dition, the resistance of the sample decreases as the car- 
rier concentration increases. In Fig. [7^ the Dirac point is 
observed to be at ^ 21 V and above this value the Fermi 
energy is shifted to the conduction band and thus the 
charge carriers are primarily electrons while below 21 V, 
the Fermi energy is shifted to the valence band and the 
charge carriers are mostly holes. Thus, this figure con- 
firms that our junctions are capable of carrying a bipolar 
supercurrent. It is also clear that the gate voltage can be 
used to tune the critical current of our graphene proxim- 
ity junctions, which can be used as a tool to study the 
escape dynamics. 

Current-phase relation in graphene and Ic{T).~ Al- 
though detailed information about CPR in the SGS prox- 
imity junction was not needed for our analysis of current 
switching statistics, nevertheless, we provide such infor- 
mation for completeness. The equilibrium Josephson cur- 
rent can be found from the thermodynamic relation 



2e dF 
h d(f> 



(7) 



by knowing free energy F as a function of the phase dif- 
ference across the junction. Bardeen et. al. !31 derived 
general expression for F with the help of Bogolubov-de 
Gennes (BdG) equations, which reads 

F = -2T J2 ln[2 cosh(£„/2T)] + f d'^r\A\^/V . (8) 

This result is a generalization of that used in BCS model 
for the case of A = const. The sum over n is just what 
would be obtained for the free energy for an assembly 
of independent fermions with energies e„. As in the 
Hartree-Fock approximation, this counts the interaction 
energy twice. The next term, J d'^r\A\'^ /V, the negative 
of the interaction energy, corrects for this double count- 
ing. Since the bulk energy oc |Ap is independent of the 
phase (p one obtains from Eq. ([8| after differentiation 



4e ■r-^ dSn ^ , e. 



2T 



8eT 



1 pOD 

- de\n 


2 cosh — 


J A 


L 2T\ 



dp 



(9) 

where we have rewritten sum over n as a sum over the 
discrete positive eigenvalues e„((?!i) (n = 1, 2, ...) of BdG 
equations, and an integration over the continuous spec- 
trum with density of states p{e, </>). The additional factor 
of 2 in the above formula, as compared to the conven- 
tional expression, accounts for two valleys in graphene. 
Although above the gap states, £ > A, do contribute to 
the Josephson current [i] their role, in fact, is subleading 
as compared to the contribution coming from the An- 
dreev bounds states with energies below the gap £„ < A. 



One thus only need to consider the first term in Eq. (|9|. 
In the context of graphene the spectrum of Andreev lev- 
els can be found from the BdG equations in the following 
form [5] 



[{ihvcr • c? — ^) ® T^Ji/i = e^/j , 



(10) 



where V''^ — (V'eiV'ft) is electron(hole) wave functions 
combined in the vector, dot product is defined in the 
conventional way cr ■ d = axdx + cFydy^ symbol stands 
for the direct product between matrices, which operate 
in the isospin (ct) and electron-hole (r) spaces. Chemical 
potential p is measured with respect to Dirac point, so 
that p = corresponds to undoped graphene. Electron- 
like and hole- like wave functions are related to each other 
at the SG interfaces r± = (±L/2,y) via specular reflec- 
tion Andreev processes. Mathematically this can be rep- 
resented as follows |5] 

Mr-) - U{e)Mr-) , Mr+) = U-\e)Mr+) , 

(11a) 

?7(£) = exp(-z<?;>/2 + i/3CTj;) , /3 = arccos(£/A) . (lib) 

Assuming hard wall boundary conditions in the y- 
direction fcy-component of particle wave vector becomes 
quantized ky = Qn = {n + 1/2)'k/W. Performing then 



Fourier transform in Eq. (10 1, matching electron/hole 



plane waves at both interfaces with the help of bound- 



ary conditions ( 11 ) and setting determinant of the corre- 



sponding matrix equation to zero one finds 

cos 2/3 



COS( 



cos^ X 



cos^ 7 



sin^Xtan^7 (12) 



which determines dispersion relation for the Andreev lev- 
els, where x — kL, k = {p,/hv) cos 7, 7 — a.rcsin{hvqn/ p). 
The last equation can be resolved analytically for e as 
the function of the channel index n and superconducting 
phase difference 4> in the form 



e„(0) - A^/l- |i„Psin2(</,/2) , (13) 

where |i„p has meaning of the transmission coefficient in 
the n*'' transversal channel, which is given explicitly by 

— 2 ' (14) 

cos(K„i) + (/i//iw)2 sin (KnL) 



with Kn — \/ (p/hv)^ — q^. Note that at the Dirac point, 
/i —7- 0, all channels are evanescent, since k„ — >■ iq^ and 
thus 

\tn\'= J. (15) 

cosh (g„L) 

Having determined the spectrum of energy states below 
the gap we can return to Eq. ([9| and find Josephson 
current in the form 16 



eA^ ^ \t,i\^sm<j) ^^^^en{(t)) 



2T 



(16) 
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where £„(</)) should be taken from Eq. ( 13 ). In the vicin- 
ity of the neutrahty point, /i <C Eth, where Exh = hv/L 
is baUistic Thouless energy, to the good approximation 
one can use Eq. (15) for the transmission coefficient. 
Furthermore, if the aspect ratio of graphene sheet is 
such that W ^ L then summation over discrete n 
can be replaced by the integration ^„ ^ dx 
with g„L — >■ X. Introducing also dimensionless variable 



z — a/ 1 — sin (0/2) / cosh x one arrives at 



2eAW 
nhL 



dz . 

cos(0/2) V' 



cos((/)/2) tanh 



2T 



COs2((/)/2) 



(17) 



There is no close analytical expression for this integral, 
except for the zero temperature limit, when tanh 
In that case [3 [7] 



2T 



^ 1. 



Is{^) 



nhL 



cos((/)/2) In 



sm 



1 - sin((/i/2) J ' 



(18) 



which coincides with the result of Kulik-Omelyanchuk for 
the case of disordered SNS junction 7j. This is rather pe- 
culiar result since calculation was done for the manifestly 
ballistic limit of graphene. From the CPR we can restore 
free energy barrier for phase slips and thus hight of the 
Gibbs potential (see Eq. 3 in the main text): 



2cl^^ AW 



eInR 



C^N 



(19) 



where we have used Ic = ci%^, \dps\4>=<l>c 
R-^ = Ae'^W/nhL with ci = 1.33, C2 = 1.08, and 
(f>c = 1-97 which corresponds to the maximum of /g (</>). 
For completeness we mention that in the longer junctions, 
L 3> W, all transmissions are exponentially suppressed. 



\tn\' 



such that summation reduced to the ge- 



ometrical progression to the leading order in with 
the result 



eA 



tanh 



(20) 



Note that for our samples the aspect ratio is at least 



W/L = 40 so that Eqs. (|T8])-([T9]) should apply near the 
neutrality point. 

At the relatively high doping when Erh ^ ^ A full 
expression for i„ is needed for the calculation of Is {<!>)■ 
Although there is no close analytical expression for the 
CPR in this case, critical Josephson current can be esti- 
mates as 



Ir 



eA 



ch 



ch 



nhv 



(21) 



where Nch has meaning of the number of propagating 
transversal channels. All above considerations apply, 
strictly speaking, to the short junctions L ^ ^ in the 



sense of the proximity effect, where ^ is superconducting 
coherence length. Note, however, that our SGS devices 
are rather at the crossover between long and short limits 
with the typical ratio L/(_ ^ 2. 

Surprisingly, we do not find solid experimental evi- 
dence in support of ballistic transport in the SGS prox- 
imity junctions. All aspects of our data, and in particular 
temperature and gate voltage dependence of the critical 
current, are in fact in good quantitative agreement with 
the predictions of diffusive SNS junction model [HIHlin]- 
Specifically we use 



cIcRn = 647rr 



E 



A\L/L,)exp{-L/L,) 



(22) 



to fit the experimental data (see Fig. 3c and 3d of the 
main text). InEq. ([22]) e„ = (2n+l)7rT, E,, = ^A^ + 
and Le = yjD/2en- In theory Eq. ([22| applies at T > 
Exh: where Eth — D/Li^ is diffusive Thouless energy 
and D — w/e/2 is diffusion coefficient. For the typical 
parameters of our samples v ^ lO^m/s, 1^ ^ 10 nm, 
L ^ 300 nm one finds Exh ^ 0.1 K which corresponds to 
T > Eth for the working temperature regime T = 0.3 K. 
Note also that as shown in the extensive study of Ref . [S] , 
Eq. (22 1 works well in the rather wide temperature range. 



At the lowest temperatures, T ^ ETh^ critical current 
saturates with the exponential accuracy [H [9] 

elcRN ^ aETh [1 - b exp{-cETh/T)] , (23) 

where a = 10.8, b ~ 1.3, and c ~ 3.4. Finally, one should 
note that product cIcRn (up to a numerical factor of 
order one) sets the magnitude of the Gibbs barrier for 
the phase slips. For example, at T ~ 0.4 K in the case of 
sample 105 taking the corresponding values Ic ^ 10 jiA 
and Rm ^ 10 one estimates Gc ^ cIcRn ^ 1 K, such 
that Gc ^ T in agreement with our earlier discussions 
in the main text. 
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